Shank3 Modulates Sleep and Expression of Circadian Transcription Factors differential expression
February 28, 2019
Description
This is the report of the analysis made for the paper Shank3 Modulates Sleep and Expression of Circadian Transcription Factors by Ashley M. Ingiosi, Taylor Wintler, Hannah Schoch, Kristan G. Singletary, Dario Righelli, Leandro G. Roser, Davide Risso, Marcos G. Frank and Lucia Peixoto.
Autism Spectrum Disorder (ASD) is the most prevalent neurodevelopmental disorder in the US that often co-presents with sleep problems. Sleep impairments in ASD predict the severity of ASD core diagnostic symptoms and have a considerable impact on the quality of life of caregivers. However, little is known about the underlying molecular mechanism(s) of sleep impairments in ASD. In this study we investigated the role of Shank3, a high confidence ASD gene candidate, in the regulation of sleep. We show that Shank3 mutant mice have problems falling asleep despite accumulating sleep pressure. Using RNA-seq we show that sleep deprivation doubles the differences in gene expression between mutants and wild types and downregulates circadian transcription factors Per3, Dec2, and Rev-erb\(\alpha\). Shank3 mutants also have trouble regulating locomotor activity in the absence of light. Overall, our study shows that Shank3 is an important modulator of sleep and circadian activity.
Differential Expression Analysis
Importing data
Importing data and filtering out those genes with cpm lesser than 1. We use the filtered.data method in NOISeq package.
countMatrix <- ReadDataFrameFromTsv(file.name.path="./data/refSEQ_countMatrix.txt")## ./data/refSEQ_countMatrix.txt read from disk!# head(countMatrix)
designMatrix <- ReadDataFrameFromTsv(file.name.path="./design/all_samples_short_names.tsv")## ./design/all_samples_short_names.tsv read from disk!# head(designMatrix)
filteredCountsProp <- filterLowCounts(counts.dataframe=countMatrix,
is.normalized=FALSE,
design.dataframe=designMatrix,
cond.col.name="gcondition",
method.type="Proportion")## features dimensions before normalization: 27179## Filtering out low count features...
## 14454 features are to be kept for differential expression analysis with filtering method 3Plot PCA of log unnormalized data
PCA Plot of filtered not-normalized data.
PlotPCAPlotlyFunction(counts.data.frame=log1p(filteredCountsProp),
design.matrix=designMatrix,
shapeColname="condition", colorColname="genotype", xPCA="PC1", yPCA="PC2",
plotly.flag=TRUE, show.plot.flag=TRUE, prefix.plot="Prop-Un-Norm")## [1] TRUEControl Genes
Negative control genes
Loading Negative Control Genes to normalize data
library(readxl)
sd.ctrls <- read_excel(path="./data/controls/Additional File 4 full list of BMC genomics SD&RS2.xlsx", sheet=1)
sd.ctrls <- sd.ctrls[order(sd.ctrls$adj.P.Val),]
sd.neg.ctrls <- sd.ctrls[sd.ctrls$adj.P.Val > 0.9, ]
sd.neg.ctrls <- sd.neg.ctrls$`MGI Symbol`
sd.neg.ctrls <- sd.neg.ctrls[-which(is.na(sd.neg.ctrls))]
int.neg.ctrls <- sd.neg.ctrls
int.neg.ctrls <- unique(int.neg.ctrls)
neg.map <- convertGenesViaMouseDb(gene.list=int.neg.ctrls, fromType="SYMBOL",
"ENTREZID")
# sum(is.na(neg.map$ENTREZID))
neg.ctrls.entrez <- as.character(neg.map$ENTREZID)
ind.ctrls <- which(rownames(filteredCountsProp) %in% neg.ctrls.entrez)
counts.neg.ctrls <- filteredCountsProp[ind.ctrls,]positive control genes
Loading Positive Control Genes to detect them during the differential expression step.
## sleep deprivation
sd.lit.pos.ctrls <- read_excel("./data/controls/SD_RS_PosControls_final.xlsx",
sheet=1)
colnames(sd.lit.pos.ctrls) <- sd.lit.pos.ctrls[1,]
sd.lit.pos.ctrls <- sd.lit.pos.ctrls[-1,]
sd.est.pos.ctrls <- read_excel("./data/controls/SD_RS_PosControls_final.xlsx",
sheet=3)
sd.pos.ctrls <- cbind(sd.est.pos.ctrls$`MGI Symbol`, "est")
sd.pos.ctrls <- rbind(sd.pos.ctrls, cbind(sd.lit.pos.ctrls$Gene, "lit"))
sd.pos.ctrls <- sd.pos.ctrls[-which(duplicated(sd.pos.ctrls[,1])),]
sd.pos.ctrls <- sd.pos.ctrls[-which(is.na(sd.pos.ctrls[,1])),]Normalizations
TMM Normalization
Normalizing data with TMM, as implemented in edgeR package, and plotting a PCA and an RLE plot of them.
normPropCountsUqua <- NormalizeData(data.to.normalize=filteredCountsProp,
norm.type="tmm",
design.matrix=designMatrix)
PlotPCAPlotlyFunction(counts.data.frame=log1p(normPropCountsUqua),
design.matrix=designMatrix, shapeColname="condition",
colorColname="genotype", xPCA="PC1", yPCA="PC2",
plotly.flag=TRUE, show.plot.flag=TRUE,
prefix.plot="TMM-Norm")## [1] TRUEpal <- RColorBrewer::brewer.pal(9, "Set1")
plotRLE(as.matrix(normPropCountsUqua), outline=FALSE, col=pal[designMatrix$gcondition])
TMM + RUVs Normalization
Applying a RUVs method of RUVSeq package on normalized data, in order to adjust the counts for the unwanted variation. And of corse we plot a PCA and an RLE plot on these data.
library(RUVSeq)
neg.ctrl.list <- rownames(counts.neg.ctrls)
groups <- makeGroups(designMatrix$gcondition)
ruvedSExprData <- RUVs(as.matrix(round(normPropCountsUqua)), cIdx=neg.ctrl.list,
scIdx=groups, k=5)
normExprData <- ruvedSExprData$normalizedCounts
ggp <- PlotPCAPlotlyFunction(counts.data.frame=log1p(normExprData),
design.matrix=designMatrix, shapeColname="condition",
colorColname="genotype", xPCA="PC1", yPCA="PC2",
plotly.flag=FALSE, show.plot.flag=FALSE, save.plot=FALSE,
prefix.plot=NULL)## [1] FALSEggplotly(ggp)dir.create("plots")
save_plot(filename="plots/PCA.pdf", plot=ggp)
pal <- RColorBrewer::brewer.pal(9, "Set1")
plotRLE(normExprData, outline=FALSE, col=pal[designMatrix$gcondition])
edgeR Differential Expression Analysis
Making differential expression analysis with edgeR package on four different contrasts.
Here is a brief legend:
- WTHC5: Wild Type Home Cage Control 5 days
- WTSD5: Wild Type Sleep Deprivation 5 days.
- KOHC5: Knock Out Home Cage Control 5 days.
- KOSD5: Knock Out Sleep Deprivation 5 days.
padj.thr <- 0.05
venn.padgj.thr <- 0.1
desMat <- cbind(designMatrix, ruvedSExprData$W)
colnames(desMat) <- c(colnames(designMatrix), colnames(ruvedSExprData$W))
cc <- c("S3HC5 - WTHC5", "S3SD5 - WTSD5", "WTSD5 - WTHC5")
rescList1 <- applyEdgeR(counts=filteredCountsProp, design.matrix=desMat,
factors.column="gcondition",
weight.columns=c("W_1", "W_2", "W_3", "W_4", "W_5"),
contrasts=cc, useIntercept=FALSE, p.threshold=1,
is.normalized=FALSE, verbose=TRUE)
names <- names(rescList1)
rescList1 <- lapply(seq_along(rescList1), function(i)
{
attachMeans(normalized.counts=normExprData, design.matrix=desMat,
factor.column="gcondition", contrast.name=names(rescList1)[i],
de.results=rescList1[[i]])
})
names(rescList1) <- namesShank3 Home Cage control VS Wild Type Home Cage Controls
volcano plot
A volcano plot of differential expressed genes.
res.o.map1 <- convertGenesViaMouseDb(gene.list=rownames(rescList1[[1]]),
fromType="ENTREZID")
res.o <- attachGeneColumnToDf(mainDf=rescList1[[1]],
genesMap=res.o.map1,
rowNamesIdentifier="ENTREZID",
mapFromIdentifier="ENTREZID",
mapToIdentifier="SYMBOL")
WriteDataFrameAsTsv(data.frame.to.save=res.o,
file.name.path=paste0(names(rescList1)[1], "_edgeR"))
vp <- luciaVolcanoPlot(res.o, prefix=names(rescList1)[1],
positive.controls.df=NULL,
threshold=padj.thr)
ggplotly(vp)de <- sum(res.o$FDR < padj.thr)
nde <- sum(res.o$FDR >= padj.thr)
detable <- cbind(de,nde)
rownames(detable) <- names(rescList1)[1]
ddetable <- detable
tot.ctrls <- dim(sd.pos.ctrls)[1]
idx.pc <- which(tolower(res.o$gene) %in% tolower(sd.pos.ctrls[,1]))
tot.pc.de <- sum(res.o$FDR[idx.pc] < padj.thr)
tot.pc.nde <- length(idx.pc) - tot.pc.de
wt <- res.o[which(res.o$FDR < padj.thr),]
wt.sign.genes.entrez <- rownames(res.o)[which(res.o$FDR < venn.padgj.thr)]
kowthc5 <- res.o[which(res.o$FDR < padj.thr),]
kowthc5.sign.genes.entrez <- rownames(res.o)[which(res.o$FDR < venn.padgj.thr)]Shank3 Sleed Deprivation VS Wild Type Sleep Deprivation
volcano plot
A volcano plot of differential expressed genes.
rs2.o.map <- convertGenesViaMouseDb(gene.list=rownames(rescList1[[2]]),
fromType="ENTREZID")
res.rs2.o <- attachGeneColumnToDf(mainDf=rescList1[[2]],
genesMap=rs2.o.map,
rowNamesIdentifier="ENTREZID",
mapFromIdentifier="ENTREZID",
mapToIdentifier="SYMBOL")
WriteDataFrameAsTsv(data.frame.to.save=res.rs2.o,
file.name.path=paste0(names(rescList1)[2], "_edgeR"))
vp <- luciaVolcanoPlot(res.rs2.o, positive.controls.df=NULL,
prefix=names(rescList1)[2],
threshold=padj.thr)
ggplotly(vp)